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ABSTRACT 

Relativistic shocks can accelerate particles, by the first order Fermi mechanism, 
which then emit synchrotron emission in the post shock gas. This process is of par- 
ticular interest in the models used for the afterglow of gamma ray bursts. In this 
paper we use recent results in the theory of particle acceleration at highly relativistic 
shocks to model the synchrotron emission in an evolving, inhomogeneous and highly 
relativistic flow. We have developed a numerical code which integrates the relativistic 
Euler equations for fluid dynamics with a general equation of state, together with a 
simple transport equation for the accelerated particles. We present tests of this code 
and, in addition, we use it to study the gamma ray burst afterglow predicted by the 
fireball model, along with the hydrodynamics of a spherically symmetric relativistic 
blastwave. 

We find that, while, broadly speaking, the behaviour of the emission is similar to 
that already predicted with semi-analytic approaches, the detailed behaviour is some- 
what different. The "breaks" in the synchrotron spectrum behave differently with time, 
and the spectrum above the final break is harder than had previously been expected. 
These effects are due to the incorporation of the geometry of the (spherical) blastwave, 
along with relativistic beaming and adiabatic cooling of the energetic particles leading 
to a mix, in the observed spectrum, between recently injected "uncoolcd" particles 
and the older "cooled" population in different parts of the evolving, inhomogeneous 
flow. 
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1 INTRODUCTION 

Relativistic shock fronts arise in astrophysics when a rel- 
ativistic flow propagates into ambient material such as in 
the jets of active galactic nuclei (AGNs) and fireball models 
of gamma ray bursts (GRBs). These shock fronts are also 
sites of energetic particle acceleration which in turn leads 
to non-thermal emission processes such as synchrotron ra- 
diation and inverse Compton scattering. Particle accelera- 
tion at relativistic shock fronts can occur via the first or- 
der Fermi mechanism. In this picture, energetic particles of 
speed V are scattered back and forward across the shock 
front by magnetic turbulence. Each crossing involves the 
particle meeting the flow head-on leading to an increase of 



tering the particle distribution is almost isotropic both up- 
stream and downstream. This results in a power law index 
for the particle spectrum which only depends on the shock 
compression ratio, r = Ui/U^- In relativistic flows the sit- 
uation is more complicated. When v is comparable to Ui,2 
the timescale for isotropisation of the energetic particle pop- 
ulation by the magnetic turbulence (which determines the 
degree of anisotropy of the particle distribution upstream 
and downstream) becomes important. Moreover even if the 
particle distribution were isotropic in one fluid frame it will 
be anisotropic when viewed in the shock or other fluid frame. 
As a result, unlike the case for nonrelativistic flows, the par- 
ticle spectrum depends not only on the shock strengt h but 



energy. 



AtT- 



uuuiclaLivisLic shucks, v >> Ui^i whcic Ui aud 
U2 are ' the upstream and downstream flow speeds in the 
shock rest frame respectively. With diffusive particle scat- 



also on the d etails of how the particles are scattered (Kirk 



Duffy 1999). This creates difficulties when studying shock 



acceleration in fiows such as spherical, relativistic blastwaves 
(Blandford & McKee 1976), the energetic particle distribu- 
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tion is the solution of a kinetic equation which describes 
scattering in the local fluid rest frame. 
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However, there has been recent progress with ultra- 
relativistic shocks, V > 10, for which it has been shown that 
the index of the power-law distribution of energetic parti- 
cles, which undergo isotropic scattering, is q = 2.23 where 
N{E) oc E^'' is the differ ential number of particles with 
energy E (Kirk et al. 2000). This result has potentially im- 



portant implications for fireball models of GRBs (Cavallo 
& Rees 1978; Goodman 1986; Paczyhski 1986) in which a 
large amount of energy is suddenly released in the form of 
an optically thick e lectron-positron plasma . The presence 
of baryonic matter ( Rees fc Meszaros 1992 ) results in the 



conversion of the initial energy into a relativistically mov- 
ing bulk flow. The forward and reverse shocks will accel- 
erate particles by the above mechanism in the presence of 
magnetic turbulence. Similarly, internal shocks can occur be- 
cause of variability in the central engine of a GRB and these 
are also sites of relativistic shock acceleration. In this paper 
we present a model for the synchrotron emission of shock 
accelerated particles downstream of the forward and reverse 
shocks in a spherically symmetric relativistic blastwave. In 
contrast to the strict fireball model our initial conditions 
are for a gas at rest with a large amount of internal ther- 
mal energy confined to a small volume. As in the fireball 
case the gas will accelerate to the point where all of the 
initial energy is converted into bulk kinetic energy if the ex- 
ternal medium is sufficiently tenuous. Ultimately the flow 
will become self-similar when the amount of swept up mass 
becomes significant and the flow will decelerate thereafter. 

Our principle aim in this paper is to study the influence 
of an evolving, inhomogeneous flow on the energetic parti- 
cles produced by the first order Fermi mechanism at rela- 
tivistic shocks. This is of central importance in the theory 
of GRBs and many papers have recently approached related 
problems on a number of different levels. The adiabatic evo- 
lution of a GRB interacting with an external medium, and 
the ass ociate d emission, has been considered (Panaitescu & 
Kumar 200C ) where an analytic model is introduced for the 



Lorentz factor of the shock front and relativistic remnant. 
A similar hydrodynamical model has also been employed in 
Moderski, Sikora & Bulik (2000) for the deceleration of the 
blast wave but for beamed ejecta. Such a geometry h as also 
been discussed in a recent paper (Granot et al. 2001) where 
a full hydrodynamical simulation is carried out. In this pa- 
per we solve the full hydrodynamical equations for a spheri- 
cally symmetric system which starts from rest but begins to 
expand relativistically, and ultimately becomes self-similar, 
as a result of the enormous amount of internal energy in 
the central engine. Moreover, particles are accelerated at all 
shock waves which arise in our simulation, which in this case 
are the forward and reverse shocks, with a spectrum deter- 
mined by acceleration theory. This model will be generalised 
to jet type geometries in future work. 

In order to achieve this we have developed a hydro- 
dynamical code, described in Sect. ^ and Appendix ^ 
which integrates the relativistic Euler equations reliably for 
Lorentz factors of up to several hundred. This code is ap- 
plicable to a gas with a general equation of state in one 
dimension and in corp orates adaptive hierarchical mesh re- 
finement. In Sect. |A1| we present tests of our hydrodynami- 
cal code. Our model for the injection of accelerated particles 
into the flow, and their subsequent cooling, is presented in 
Sect, pi The initial conditions are presented in Sect. W, while 



in Sect. ^ we give a discussion of the hydrodynamics of blast- 
waves with non-zero initial radius. The results are presented 
and discussed in Sect. 0. 



2 HYDRODYNAMICS 

We wish to simulate a spherically symmetric explosion which 
gives rise to relativistic velocities. We use the relativistic Eu- 
ler system of equations in order to calculate the evolution of 
this system. The code used is a second-order Godunov-type 
scheme (e.g. van Leer 1977) which uses a linear Riemann 
solver for shocks, and a non-linear one for strong rarefac- 
tions. Appendix N gives the details of the code, along with 
some test calculations. In this section we briefly describe the 
equations used. 

The conservation equations for inviscid relativistic hy- 
drodynamics in spherical symmetry are 
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where F is the fluid Lorentz factor, p is the proper density, 
P is the velocity in c = 1 units, w is the enthalpy and p is 
the proper pressure. Time, t, and distance, r, refer to the 
coordinates measured in the observer's frame. We can relate 
the enthalpy, density and pressure by 

7* 



w = p~\- 



7. 



(4) 



If the gas is composed of A'^ species which each ha ve mass 
rrii and number density Ui then, for a Synge gas (Falle & 
Komissarov 1996), 



7* 



^2 (CO 
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(5) 



where = and Ka are the modified Bessel functions. 
In this work we assume two species - electrons and protons 
- which are in thermal equilibrium. The variable 7* defined 
above is different to the ratio of specific heats, 7 which is 
used in the calculation of, for example, sound-speeds. With 
this relation we should be able to solve the evolution equa- 
tions |l| to ^ However, it should be noted that in relativis- 
tic hydrodynamics there are no explicit relations giving the 
primitive variables p, j3 and p in terms of the conserved 
variables Fp, wV^ (3 and wV^ — p. Therefore an iterative al- 
gorithm must be developed to do this. The precise nature of 
this algorithm can affect how high a Lorentz factor can be 
simulated by the code. We have found that the code used 
here will simulate flows up to Lorentz factors of at least 
several hundred reliably. 



3 PARTICLE ACCELERATION AND 
SYNCHROTRON EMISSION 

3.1 The particle distribution 

Relativistic shocks, in the presence of magnetic fluctuations 
which enable multiple shock crossings, can accelerate par- 
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tide s by the diffusive shock mechanism (Kirk & Schneider 
198S |). -| Vs described in the introduction it has been shown 
that at ultra-relativistic shocks this produces a power law 
distribution of particles with index 2.23, i.e. N{E) oc E~^-^^ , 
so that the energy density does not diverge at high particle 
energies. The acceleration timescale in relativistic flows will 
scale roughly with the particle mean free path in the tur- 
bulent magnetic field divided by the shock speed; although 
there is no published derivation of the precise acceleration 
timescale. In this paper we are interested in what happens 
to the power law population produced by highly relativistic 
shocks over the much longer timescale which is the size of 
the blastwave divided by c. Therefore, over hydrodynamical 
timescales, we can treat acceleration as an impulsive injec- 
tion of energetic particles, with this universal power law up 
to arbitrarily high energies, into each fluid element which 
passes through a relativistic shock. We will address the role 
of an intrinsic cut-off in the shock accelerated population, 
either as a result of a finite acceleration timescale or loss pro- 
cesses, in a future paper. The particles subsequently suffer 
synchrotron and adiabatic losses from which we can calcu- 
late the emission. We need to introduce three free param- 
eters which, when taken together with the results of our 
hydrodynamical simulations, will give a complete picture of 
the blastwave's evolution and the associated emission of en- 
ergetic particles. These parameters are 

• the ratio, e;,, between the magnetic field energy density 
and the thermal energy density, 

• the fraction, te, of the downstream thermal energy den- 
sity which is converted into high energy electrons as a fiuid 
cell is shocked and, 

• the energy, -Bmin, of the lowest energy particle which is 
produced at the shock. The differential energy spectrum is 
a power law in energy of index 2.23 from i?min. 

Once an electron is injected at the shock it will be scat- 
tered in the local fiuid frame and suffer both synchrotron and 
adiabatic losses. If the length scale over which the electrons 
are isotropised is much shorter than any other length scale of 
interest then the relativistic electrons will respond adiabat- 
ically to the expansion or contraction of the fiow. Adiabatic 
losses, or indeed gains, are then described by the fact that 
p/p^^^ is constant where p and p are the particle momen- 
tum and fiuid density in the local fiuid frame. With E — pc 
for ultra-relativistic particles the combined synchrotron and 
adiabatic losses are described in the comoving frame by 



E = -aB'^E'^ + l:-E 



3p 



(6) 



where a is a constant and B the magnetic field strength in 
the local fiuid frame. Consider now a fiuid element which is 
shocked at time to when a power law distribution of energetic 
particles is injected and where t is time measured in the 
comoving frame. This population then evolves according to 
the equation 



(7) 



where N{E, t) is the differential number of particles of en- 
ergy E at time i. The losses, E, are given by equation 
while the injection term is Q{E,t) — QoS{t — to)H{E — 
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EminjE"^ [H being the Heviside function) which describes 
injection of a power law spectrum, starting at time to with a 
minimum particle energy i^min. We can solve equation (Q) by 
finding the characteristic curves along which N{E,t)dE — 
N{Eo,to)dEo where a particle with energy Eq at to cools to 
an energy E at time i. From equation M we have 



d_ 

di 
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which can be solved to show that along a characteristic curve 



Eq 

p^ ~aE p^B'^di 
N is conserved so that 
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and the solution to (0) becomes 



N{E,i) =QoPo^ p^E-^[p^-aE / p^B'^dt 



p-2 



(11) 



In order to calculate the particle spectrum we need to 
evaluate the integral in equation ( pT| ) for each fiuid element 
which has passed through a shock at some point. No fiuid el- 
ement passes through two or more shocks during these simu- 
lations; our method would have to be adapted in such a case 
(e.g. for overtaking internal shocks). We introduce a conti- 
nuity equation for the integral quantity in equation (pT[). 

Define / = f El^—dt where the factor of F^^ boosts ob- 
server time to the comoving fiuid time. Then the equation 
for the distribution of / is simply 



dr 



(12) 



Thus the quantity / increases by p^^^ B^ At/V with each 
time-step At, following the fiuid element, as required. If the 
gradient of F/3 is ever less than a certain negative value, cho- 
sen to ensure that the fiuid element has passed through a 
shock, then I is set to at that point. The above conserva- 
tion equation can easily be solved in conjunction with the 
hydrodynamical equations using the same numerical meth- 
ods. Note that / is not kept at zero in unshocked fiuid, but 
this does not affect the synchrotron emission since, in such 
fluid, there are no energetic particles. 



3.2 Synchrotron emission 

We approximate the emission of a single particle to be a 
delta function in frequency with a single particle emissivity 
in the fluid frame given by joi'j) = ao^^B^5{v — ai^^B) 
where 7 is the Lorentz factor of the electron in the fluid 
frame, oo = arc/Gn and ai = e/2-Rm^c with ar the Thom- 
son cross section. With the local electron spectrum given by 
[ll| , the emissivity in the local fluid frame is then 
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simple case to the treatment used in Moderski et al. (2000). 
In particular from those particles which have not radiated 
away a significant fraction of their energy we get the un- 
cooled spectrum of oc jy^tP"^)/^ while the higher energy 
particles have a steeper, cooled spectrum of Fi, oc u^^'''^. 
Below the minimum observed frequency the flux scales as 
Fi, oc ly^^^. At even lower frequencies the effect of syn- 
chrotron self-absorption will become important, although we 
have not explicitly calculated this, leading to a part of 
the spectrum. 



4 INITIAL CONDITIONS 



Figure 1. Diagram of the technique used to calculate the spec- 
trum of photons arriving at an observer at a given time t{obs). 
The fluid element is located at r in the simulation, and the cur- 
rent time is t. This, in conjunction with L, the distance between 
the observer and the blastwave, defines 6 and <f>. See text. 



(13) 



The emissivity in the observer's frame is related to that 
in the fluid frame by £^ = D^Ep. The Doppler factor, D, for 
the fluid element with a three velocity /3 making an angle 9 
to the line of sight is D = [r(l — where /i — cos 9. The 
observed frequency is related to the photon frequency in the 
fluid frame hj v = DO. The observed intensity of radiation 
at an inclination (j) to the line between the observer and the 
centre of the explosion is 



£u{s, <f), t) ds 



(14) 



where to is the time of observation and t = to — s/c is the 
time of emission in the observer's frame a distance s away 
from the observer (figure The observed flux density is 
obtained by integrating the intensity over the entire source 

F^{to)= / I.i(t>,to)dn (15) 
Jn 

In cylindrical coordinates {z,r) centred on the explosion, 
s ^ L ~ z since rA ^ 1 so that 



FAto) 



2% 

Z2 



dz / £u{z,r,t)rdr 



(16) 



with L the distance from the observer to the centre of the 
explosion and rc the size of the computational domain. If ro 
is the initial radius of the explosion then the interval of to 
since the arrival of the first photon is to = to — [L — ro)/c 
so that the flux can be written as F^{to). By the change 
of variable z = ro + c{t — to) the integration over z can be 
converted to one over t, the time of emission in the observer's 
frame, 

_ 2-Kc f'""^ r' — 

F,.{to) = jYj '^^ J £Aro + c{t~-to),r,t)rdr (17) 

In the simple case where the fluid is assumed to move 
directly towards the observer with a constant velocity we 
have successfully compared our method with the exact re- 
sults. The above anal ysis, w hich is essentially the same as in 
Komissarov & Falle ( 1997), also corresponds for the latter 



In studies of the fireball model of gamma-ray bursts, the pa- 
rameters used are the total mass, M, and the initial radius, 
-Ro- In general, E, the total energy, is thought to lie some- 
where in the region lO^^-lO^'' ergs. The ratio between the 
rest-mass energy in the blast (i.e. the mass of the baryonic 
component), and the total energy determines two critical 
radii. These are the radius at which the baryons have been 
accelerated up to their maximum velocity, i?c, and the ra- 
dius at which the shell of ejected baryons have swept up their 
own mass in interstellar material, R^. This latter radius is 
the radius at which the baryonic shell begins to decelerate. 
The former also depends on _Ro, the initial radius of the 
blast. This initial radius is thought to be quite small, with 
Ro ~ 10^^ cm. 

In order to get all the initial energy in the blast con- 
verted into kinetic energy of the baryons, it is also required 
that the fireball be optically thick to pair creation until the 
radius of the blast exceeds Rc- If this were not the case, then 
the photons would escape into space before accelerating the 
baryons to their maximum velocity. In our case we put the 
initial energy of the blast into thermal pressure. This means 
that it is already in the kinetic energy of the baryons. Hence 
the latter consideration need not be taken into account in 
our initial conditions. 

The results of two simulations are presented in this pa- 
per. Each of the simulations used the following properties: 



580. This value 



• E^ 10^^ ergs 

• Ratio of energy to mass: r] = 
being chosen so that the Newtonian value of Rd is twice 
Rc- In reality, Rd is slightly less than R^ due to relativistic 
effects. 

• Ro ~ 1.2 X 10^* cm. This is much larger than the postu- 
lated physical value. However, reducing Ro to 10^^ cm would 
impose extremely severe computational costs on the simu- 
lation. This increase in the value of _Ro is not expected to 
change the properties of the synchrotron emission calculated 
as we are interested in the afterglow in this work, and not 
the emission from the acceleration phase. 

• £e = 0.01 - this value is consistent with observed values 



found by, e.g. Wijers & Galama (199£) 

• ei, = 0.01 in one simulation, and e;, = 0.1 in the other. 
This is the only way in which the two simulations differ. 

We must also choose the injection energy -Emin for the 
shock acceleration process. Physically this will be set by non- 
linear plasma processes. The highest energy thermal parti- 
cles downstream of a shock can leak into the upstream form- 
ing a beam of particles which can excite plasma turbulence. 
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This turbulence then scatters particles back into the down- 
stream fluid and multiple shock crossings become possible. 
The problem of determining the details of this process is an 
unsolved one in the theory of particle acceleration at rela- 
tivistic shocks. However, it is clear that the value Emin must 
correspond to the high energy tail of the downstream ther- 
mal population. This is hardly surprising since the thickness 
of a shock wave will be set by the gyroradius of a thermal 
particle in a coUisionless plasma and first order Fermi ac- 
celeration will only be applicable to particles with gyroradii 
much larger than this. Consequently we have fixed the value 
of Emin = IQmeC? which is chosen so that it exceeds the 
thermal energy for all shocked fiuid elements. Setting -Emin 
as a factor of a few times the thermal energy at the instant 
a piece of fiuid is shocked does not change our results ap- 
preciably. 

Initially, then, there is a sphere of radius Rq which is at 
rest at the origin. This sphere contains all the initial energy 
and mass of the blast. Outside this sphere is a gas with a 
density of 1 cm""^. 



5 HYDRODYNAMIC EVOLUTION OF A 
SPHERICAL BLASTWAVE 

In this section we discuss the qualitative evolution of a spher- 
ical blastwave with finite initial radius. This discussion ap- 
plies to both relativistic and non-relativistic cases. The pur- 
pose is to elucidate, in a qualitative fashion, how the for- 
ward and the reverse shocks are formed, and how they be- 
have, from a hydrodynamic perspective. Figure ^ shows the 
early-time evolution of the system, computed using the code 
presented here. The initial conditions used were as follows: 
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Figure 2. Plots of the spatial component of the 4-velocity, along 
with the proper density for times of t = 1.5 (top) and t = 3 
(bottom). All the waves in the system are labelled. See text. 
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on r < 1 
on r > 1 



These initial conditions are similar to those in Sect. 



A1.3, and are chosen so that the relevant features are clearly 
visible in the results. 

Initially a shock is driven into the ambient medium (the 
"forward shock" in figure |^ , while a rarefaction propagates 
from r = Ro back to r = (the "ingoing rarefaction" in 
figure ^). Note that, at t — 1.5, the ingoing rarefaction has 
still not reached the origin. The forward shock gets progres- 
sively stronger as blast material gets accelerated down the 
pressure gradient of the rarefaction. There is also an entropy 
wave (the "contact discontinuity") which moves out from Ro 
with a speed slightly lower than that of the blastwave. 

In addition, a reverse shock is created almost imme- 
diately between the entropy wave and the rarefaction (see 
figure If the system were in planar symmetry then one 
would not expect this reverse shock since there are only 
three characteristics in hydrodynamics (see, e.g. Landau & 
Lifschitz 1966). We get this fourth wave purely as a result 
of the non-planar nature of the system. It forms because, as 
the forward shock moves outwards it sweeps up an increas- 
ing amount of interstellar matter per unit distance. This 
means that the forward shock does not move outwards as 



fast as it would in the planar case, and the reverse shock is 
formed to decelerate material behind the forward shock to 
the appropriate speed. 

When the rarefaction reaches r = it is reflected and 
enhanced. This "reflected rarefaction" can be seen clearly at 
t = 3 in figure This rarefaction will eventually catch up 
with the reverse shock/blastwave system (see figure ^ and 
Sect. |6.l| ). When the rarefaction catches up with the reverse 
shock, the ram-pressure of the material entering this shock is 
drastically reduced, because it has been rarefied and decel- 
erated. Hence it can no longer support the thermal pressure 
between the forward and reverse shocks. As a result, the 
material between these two shocks expands. Since the ram- 
pressure of the material entering the forward shock has not 
changed (assuming constant external density), this neces- 
sitates the reverse shock slowing down with respect to the 
rest-frame of the initial blast. Indeed, it will actually begin 
propagating towards the origin. The speed at which it does 
this depends on the initial value of rj, since this determines 
the strength of the initial (and so the refiected) rarefaction. 
For example, if rj is very small, then the force due to the 
pressure gradient will produce a relatively small increase in 
velocity of the material in the initial sphere. Hence the rar- 
efaction will be weak. If, on the other hand, rj is large then 
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the force due to the pressure gradient will accelerate the 
blast material to high velocities very quickly, producing a 
very strong rarefaction. For the conditions in Sect. ^ it be- 
comes sufficiently strong to cause ejected material to flow 
back towards the origin also. 

When the reverse shock reaches r = it rebounds and, 
as it propagates outwards again, it weakens. Once this shock 
is sufficiently weak to be ignored, we have reached the self- 
similar stage of evolution. Prior to this point the finite ra- 
dius of the initial blast plays a part, and hence there is a 
significant length-scale in the system. If the blast was ini- 
tially of zero radius then the whole process described above 
would happen infinitely quickly (or, equivalently, would not 
happen at all). 

We can relate the above "hydrodynamical" picture with 
the one used by, for example, Rees & Meszaros (1992), as 
follows. The coasting radius, Rc, is reached when all ma- 
terial has been accelerated down the pressure gradient of 
the ingoing rarefaction. This, of course, never occurs, but 
one can define Rc as being the radius when "most" of the 
material has been accelerated. 

The deceleration radius, Rd, is then the radius at which 
the reflected rarefaction catches up with the reverse shock. 
This is when the material between the forward and reverse 
shock expands back towards the origin, necessarily deceler- 
ating significantly as it does so. 



6 RESULTS 

We discuss the results of the simulations in two sections. The 
first deals with the detailed hydrodynamic evolution of the 
blastwave. The second concerns the observed spectra and 
light-curves. 



6.1 Hydrodynamic results 

Here we show the hydrodynamic aspect of the results of our 
simulations. Figure ^ shows the distribution of proper den- 
sity, velocity and pressure after 1 x 10"" seconds. We can see 
the two rarefactions mentioned in Sect. |^. The head of the 
reflected rarefaction lies at the location of the peak velocity, 
with its tail at r = 0. The original rarefaction lies between 
the peak of velocity and the rise in the density and pressure 
plots at ~ 1 X 10^ light-seconds. 

The forward shock can be identified as the right-most 
rise in pressure and density, while the reverse shock is just 
to the left of this. 

Figure ^ shows the same plots as figure ^, but for a 
time of 7 X 10^ seconds. Here we can see the reverse shock 
has begun to separate significantly from the forward shock. 
The reflected rarefaction is now the only rarefaction present 
in the system. The reverse shock is expanding back into 
this rarefaction due to the insufficient ram-pressure of the 
material in the reflected rarefaction to support the pressure 
of the material between the forward and reverse shocks. 

It is possible to see oscillations in the density plot just 
behind the reverse shock. These oscillations arise out of nu- 
merical errors during the time that the reflected rarefaction 
comes close to, but is not in contact with, the reverse shock. 
The rarefaction is extremely strong here, and it is not sur- 
prising that numerical errors become significant. However, 
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Figure 3. Plots of density, 4-velocity and pressure (top to bot- 
tom) for a time of 1 X 10^ seconds after the initial blast. We can 
see the original and reflected rarefaction, as well as the forward 
and reverse shocks. See text. 



the influence of these errors on the overall solution, or on 
the calculated synchrotron emission, is negligible due to the 
low density of the region containing the errors. 

Figure |^ shows the latter stages of the non-self-similar 
evolution of the system. At this stage the reverse shock has 
propagated all the way to r = and has rebounded. We can 
see that it has weakened greatly, and will continue to do so. 
Once this shock is sufficiently weak to be ignored the sys- 
tem will evolve in a self-similar way. We can see the entropy 
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Figure 4. As for figure jsj but for a time of 7 X 10® seconds after 
the initial blast. 
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Figure 5. As for figure || but for a time of 1.4 X 10^ seconds 
after the initial blast. 



errors generated by the very strong reflected rarefaction re- 
maining just ahead of the reverse shock now. These errors, 
although apparently signiflcant here, do not affect the syn- 
chrotron results due to the low densities in this region, and 
the low relativistic boosting suffered by emission from this 
material. Both these effects make the influence of emission 
from this material negligible. 



6.2 Spectra and light curves 

In this section we present the results from the calculations 
of the synchrotron emission. 



6.2.1 Spectra 

Figure ^ shows a plot of the spectrum observed 24 hours 
after the initial blast could have been observed for the case 
where et = 0.1 (hereafter referred to as the high ei, case) and 
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Figure 6. The spectrum observed at 24 hrs after the initial blast 
observation for the case of t^, = 0.01 ("Low B-fields") and ej, = 0.1 
("High B-fields"). 



for £(, = 0.01 (hereafter referred to as the low e^, case). It is 
clear that in both cases the spectrum is a broken power-law. 

In the first section of the spectra, the flux goes as 
v^^^ , as expected for synchrotron radiation below the peak 
frequency of the lowest energy electron. At lower frequen- 
cies, which are not plotted here, synchrotron self-absorption 
would give a flux of u^^^ . The second part of the spectrum 
is, in both cases, the power law y''^-^^^ indicating that this 
part of the spectrum is dominated by emission from elec- 
trons, with a distribution of B"^'^^, which have not suffered 
significant adiabatic or synchrotron cooling. 

The break from this part of the spectrum to the fi- 
nal, steeper, part occurs in different places in the spectra 
plotted. For the high tb case, the break occurs at a lower 
frequency than for the low eb case as would be expected 
since, in higher magnetic fields, the losses suffered by the 
energetic population are correspondingly higher. However, 
while in the flnal part of the spectrum the exponent pre- 
dicted from simple theory would be —p/2 — —1.115, we find 
the spectrum to be slightly harder. In the high et case, we 
have Fi, oc v'^"^^, and in the low et case, oc z/""'^*^. 
The reason for the slightly harder spectra at high frequen- 
cies is the non-uniform velocity distribution in the "shell" of 
ejected material (see, e.g., figure^. Material moving at high 
velocity towards the observer will be more heavily weighted 
in the spectra than material moving at lower velocity, due 
to relativistic beaming. Since such high velocity material 
occurs immediately behind the forward shock, this material 
will contribute more to the spectrum than material further 
back from the shock. The fluid just behind the shock con- 
tains an electron population which has only recently been 
accelerated and so emission from here will be dominated by 
uncooled electrons to a very high frequency. Material from 
further behind the shock will have emission dominated by 
cooled electrons down to a lower frequency. Hence, above a 
certain cut-off, we expect to get a mixture between emis- 
sion from cooled and un-cooled electrons, with un-cooled 
electrons being preferentially weighted. This leads to an ex- 
ponent for the spectrum lying between the uncooled value 



case gives a slightly steeper spectrum at these frequencies 
than the low et case. For high magnetic fields the contri- 
bution from cooled electrons will be stronger closer to the 
shock than in the low magnetic field case. Therefore, while 
we still expect the uncooled electrons to be preferentially 
boosted in the spectrum, the effect of the cooled electrons 
will be stronger for higher magnetic fields. 

It is worth discussing the behaviour of the 'critical' fre- 
quencies in the spectra with time - i.e. where the breaks 
occur. The lower break, which arises due to the presence of 
electrons injected with the minimum energy JSmin behaves 



roughly as t 



This is the same for both the low and 



high £(, cases. This is a much less dramat ic de crease with 
time than that predicted by Sari & Piran ( 1997 ) and is due 
to the fact that, here, we define the minimum energy of injec- 
tion a priori. This energy is chosen so that the gyroradius of 
an electron at this energy would be accelerated by the Fermi 
mechanism. The latter authors choose to define the number 
of electrons accelerated, and this, in conjunction with the 
energy in energetic particles, defines Emin- This produces 
the different behaviours of the lower break in the spectra 
with time. 

We find that the frequency of the second (upper) break 
goes like t~^'^ for the low case and for the high et 

case. This is a much fas ter de crease with time than that pre- 
dicted in Sari & Piran (1997). This difference cannot be due 
to the different choice of definition of Emin, but it could be 
hypothesised that it is due to the effects of adiabatic cool- 
ing. This cooling will cause the second critical frequency to 
decrease faster where the ffow is divergent. However, simu- 
lations of this system without adiabatic cooling show that, 
while adiabatic cooling does have this effect on the observed 
spectrum, it is too small to explain the discrepancy. It is 
not clear what specifically causes t his diff erence between our 
results and those of Sari & Piran (1997), but it must be re- 



(i^ = —0.615) and the cooled value (- 



-1.115). This 



conclusion is further supported by the fact that the high et 



membered that here we perform a much fuller treatment of 
the dynamics and so it is not surprising that we get some 
significant differences with semi-analytic approaches. 



6.2.2 Light curves 

In this section we discuss the light curves resulting from the 
simulations described above. Figure Qshows plots of the light 
curves at low, medium and high frequencies. These frequen- 
cies are chosen so that the points are always in the regime 
where Fi, oc v^^^ (low frequencies), oc (medium fre- 

quencies) and, for high frequencies, the frequency is always 
in the third section of the spectrum. 

It should be noted in the following discussion that the 
range of times covered by these light curves is 1 to 24 hrs 
after the first signal of the blast reaches the observer. Ini- 
tially we see an increase in all the light curves shown, with 
the exception of the high frequency, high et case, where the 
emission is roughly constant over the first couple of hours. 
This increase is much faster than that predicted by previous 
analytic work. For example, for the low et case at low fre- 
quencies, the flux goes as roughly * at early times (cf Sari 
& Piran, 1997, where the maximum increase is t^^^). This 
effect is not due to any problems with numerical resolution 
since the emission in the 1st hour is dominated by mate- 
rial between the forward and reverse shocks approximately 



© 2010 RAS, MNRAS 000, 



Relativistic blastwaves and synchrotron emission 9 



le+09 



■g le+08 



le+07 



High B-fields 
Low B-fields 



1000 



10000 
Time (s) 



100000 



le+09 



le+08 



E le+07 



le+06 





High B-fields 
Low B-fields 























1000 



10000 
Time (s) 



100000 



1000 



100 



10 



High B-fields 
Low B-fields 




1000 



10000 
Time (s) 



100000 



Figure 7. Light curves calculated at low, intermediate, and high 
frequencies. See text. 



6 X 10® seconds after the initial blast, when the two shocks 
are well- resolved and separated by the code. 

The behaviour of the low frequency light curves is qual- 
itatively different to the other two. It keeps increasing for 
the duration of the simulation. This is unsurprising as it 
will only begin to decrease when the critical frequency Vmin 
corresponding to Fmin passes through the frequency of the 
light curve, and we have chosen this frequency so that this 
does not happen. It is very difficult to say that the light 
curve behaves anything like a power-law with breaks over 
this time-scale based on our results. 



The medium frequency light curves initially increase, 
as already stated, and then fall off dramatically. Again, it is 
difficult to see any point at which the behaviour of the light 
curve is a true power-law. The high frequency light curves 
also begin to fall off very steeply, and again, it is difficult to 
see any sign of a power-law behaviour in the curves. 

The lack of a clear power-law behaviour in the light 
curves may well be due to the restricted time-scale over 
which the curves are calculated (from 1 hour to 1 day) . How- 
ever, it is clear that, if there is a broken power-law behaviour 
then the breaks are smeared out, and, in addition, the final 
fall-ofi' of the light-curve would seem to be much faster than 
previously predicted. 

The unexpected behaviour of the low and medium fre- 
quency light curves is more likely to do with the different 
treatment of the minimum energy in the accelerated electron 
population. In our case, since we fix iJmin a priori, we allow 
the available energy to define the number of energetic elec- 
trons. Previously it has been common to define the number 
of energetic electrons and hence the available energy defines 
i^min- This will lead to quite different behaviours of the lower 
break frequency with time. 



7 CONCLUSIONS 

We have presented a model for the synchrotron emission 
of energetic particles downstream of relativistic, spherical 
shock waves. The hydrodynamical part of the problem has 
been solved numerically with the simpler simulations agree- 
ing with results published elsewhere. On the other hand the 
particle acceleration aspect has been treated as an injection 
process with relativistic shocks leaving a population of ener- 
getic particles immediately downstream. The spectral index 
is known from semi-analytic work so that we need only spec- 
ify the fraction of the downstream thermal energy which is 
converted into energetic particles and a lower cut-off en- 
ergy. The particles subsequently lose energy by synchrotron 
cooling and adiabatic losses. The hydrodynamical results 
have captured the evolution of both the forward and reverse 
shocks, which are of principal interest for particle accelera- 
tion and radiative emission, as well as the rarefaction waves. 
The spectra emitted from our system largely agree with the 
simple predictions for the low energy, uncooled part of the 
spectrum. However, at the higher frequency, cooled, part of 
the electron population the relativistic boosting of the mate- 
rial coming straight towards the observer hardens the spec- 
trum from the pure cooled value which comes from material 
further downstream and from material which is not moving 
directly towards the observer. Further, non-trivial behaviour 
has been found for the temporal variation of both the break 
frequencies and the light curves. Adiabatic losses coupled 
with integrating the spectrum over a spherical system tend 
to smear out the break frequencies which are predicted from 
simple scaling arguments. 

It is clear that we need to include several other effects 
before bringing our results into contact with observations. 
A generalisation to a 2-D hydrodynamical code has already 
begun. The role of internal shocks will also be included along 
with their effect on acceleration. The energetic particles obey 
a transport equation, the solution of which determines the 
spectral index and we are working towards including such a 
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solution in our model. More complicated radiative processes 
can also be included such as the inverse-Compton and syn- 
chrotron self-Compton processes. However, what we have 
done in this paper is to combine a realistic hydrodynami- 
cal model with results from particle acceleration theory in a 
relativistic flow. The computed hydrodynamical evolution, 
spectra and light curves will allow us to compare this model 
with observations in future work. 
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APPENDIX A: NUMERICAL DETAILS 

In this appendix we discuss the details of the code used for 
the calculations presented in this paper. We also present the 
results of three test calculations. 

We employ a second order finite volume Godunov-type 
scheme to solve the equations |l| to ^ in spherical symme- 
try. Assuming that the grid cells are defined so that cell i 
occupies the space [ri-\/2,'i'i+i/2\ then the scheme can be 
written as 



n + l/2 



3Af 



H/2 



-1/2 



2 p,n+l/2 



l/2-f^i-l/2 



•n + l/2 



At 



(Al) 



where superscripts refer to the time index and subscripts 
refer to the spatial index. Also, 



±1/2 



i±l/2 



,0 



(A2) 
(A3) 

(A4) 



It should be noted that the so urce term S should be vol ume- 
averaged over [ri_i/2, ri_|_i/2] (Falle & Komissarov 1996). So, 



for spherical symmetry, and using a first order scheme we 
get that the source term for the radial momentum equation 
is 

ri+i/2 + ri-i/2 



C-n _ Q n 



. — (A5) 

l/2+^«+l/2'"i-l/2+?'tl/2 

For a second order scheme, which allows p to vary linearly 
in a given cell the correct form for 5*" is 

c.n _ o 1 n+l/2 -I- n-i/2 

'^r+1/2 + r^+l/2r^-l/2 + ^f.^/j 



+gp^ 



9(ri+l/2 + n-l/2)^('-?+i/2 + ^^-1/2) 



(A6) 



4('-?+l/2 + r^+l/2r^-l/2 + '^Li/2)' 

where gp" is the gradient of the pressure in cell i at time-step 
n. 

The fluxes F are calculated from the primitive variables 
extrapolated to the cell edges using non-linear averaging of 
the gradients (e.g. van Leer 1977; Downes & Ray 1998), mak- 
ing the code second order in space. A non-linear Riemann 
solver is used in the presence of strong rarefactions, and 
otherwise a linear Riemann solver is used in order to save 
computational time. The flux terms F^^^j^^ are calculated 
using a first order Godunov-scheme from the conditions at 
time n. This makes the code second order in time. 

The details of the linear and nonlinear Riemann solvers 



follow closely the discussions of Falle (1991) and Falle & 
Komissarov (1996). 



The code also contains an option for hierarchical grid 
refinement. The a lgorit hm used for this is loosely based on 
that of Khokhlov ( 1998 ). However, this option is not used in 
the work presented here, so we will not discuss this further. 



Al Tests of the hydrodynamical code 

The code used in this work has been rigorously tested 
against, in particular, results presented in Falle & Komis- 
sarov (1996). The code reproduces their published results. 
Below we show three examples of these tests for complete- 



Al.l Test 1: Shock-tube problem with 7 = | 

In this test we use the following initial conditions: 



= 1 
= 



10'' 
10" 



on a:: < 0.5 
on a:: > 0.5 
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with gradient zero boundary conditions on x = and x = 1. 
This calculation was done for a uniform grid with 400 cells. 
The results are shown in figure along with the exact 
solution. The plots are of the proper density, and are at 
different times. It can be seen that, while the contact dis- 
continuity becomes quite smeared (as expected, since it is a 
linear wave) the shock and rarefaction are treated relatively 
well. 

In particular, the position of both the shock and the 
rarefaction are very well calculated, indicating that both the 
linear and non-linear Riemann solvers are working well. 

A1.2 Test 2: Shock-tube problem for a Synge gas 
In this test we use the following initial conditions: 



P 

u 



10^ 



on a; < 0.5 
on a; > 0.5 



with gradient zero boundary conditions on a; = and a; = 1. 
This calculation was done for a uniform grid with 400 cells. 
The results are shown in figure A2, along with the exact 
solution. Again, it can be seen that the numerical results 
match the analytic results reasonably well, although numer- 
ical dissipation does affect the contact discontinuity quite 
significantly. 

A 1.3 Test 3: Spherically symmetric expanding wave 

This test is very similar to the simulations used in this work. 
Therefore it is an interesting test to run. The conditions are 
as follows: 



P 

u 

P 



10 




10'' 

3 X 10" 



on r < 0.5 
on a; > 0.5 



The grid spans r = to r = 4 with 100 uniformly spaced 
cells. The results obtained are very similar to those given in 
Falle & Komissarov (1996) and are shown in figure A3. 

Overall, then, the code reproduces published results 
very well. One additional test was performed. This was a 
test of the code for the non-relativistic Sedov blastwave. The 
code successfully reproduced the behaviour of the radius of 
the blastwave (with r ex tt). 



It is noted in Sect. 3.1 



that the simulations performed 
here give entropy errors in some parts of the fiuid. These 
errors arise when the refiected rarefaction (see Sect. |^ has 
almost caught up with the reverse shock. In this case there 
is an extremely strong rarefaction, with differences of the 
Lorentz factor of the fluid of 30 - 40 between adjacent grid 
cells, just behind the reverse shock. These very extreme con- 
ditions result in negative pressures being produced (in the 
grid cell, not in the Riemann solver), subsequently resulting 
in entropy errors being apparent in the solution. 

As stated in Sect. 3.1, however, these errors occur in 



a region of the flow where the density of emitting particles 
is relatively low, and the errors do not affect the simulated 
emission appreciably. 
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Figure Al. Plots of the proper density at t =0.12, 0.24 and 
0.36. The points are the numerical solution, while the line is the 
exact solution. The rarefaction is tracked extremely well, while 
the shock is slightly smeared. The contact is more smeared, as 
expected. See text. 
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Figure A2. Plots of density, spatial component of the 4- velocity 
and pressure (top to bottom) for the shock-tube test. The points 
are the numerical solution, while the line is the exact solution. 
The numerical soluti on is v irtually indistinguishable from that in 
Falle & Komissarov (Fooi). 



Figure A3. Plots of density, spatial component of the 4- velocity 
and pressure (top to bottom) . The points show our solution, while 
the lin e shows the solution from the code of Falle & Komissarov 
(1996). It can be seen that two compare very well. 
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